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Abstract We study the linear theory of Kelvin-Helmholtz instability in a 
layer of ions and neutrals with finite thickness. In the short wavelength limit 
the thickness of the layer has a negligible effect on the growing modes. However, 
perturbations with wavelength comparable to layer's thickness are significantly 
affected by the thickness of the layer. We show that the thickness of the layer 
has a stabilizing effect on the two dominant growing modes. Transition between 
the modes not only depends on the magnetic strength, but also on the thickness 
of the layer. 
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1 Introduction 

When two fluids are in relative motion on either side of a common boundary, 
the Kelvin-Helmholtz (KH) instability can occur. This instability, the sim- 
plest example of shear flow instability, is a well-known phenomenon in fluid 
mechanics and astrophysics (e.g. Chandrasekhar 1961). For example, the in- 
teractions between the solar wind and planetary magnetospheres are studied 
based on the KH instability (e.g. Nagano 1979). The KH instability, because 
of the velocity shear in the mass flow between coronal plume structure and 
the interplume region of the sun, could be a potential source for some of the 
Alfvenic fluctuations observed in the solar wind (Andries, Tirry & Goossens 
2000). The structure and the helical wave motion in the ionized cometary tails 
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as studied by Ershkovich & Mendis (1986), Ershkovich, Prialnik & Eviatar 
(1986) and Ershkovich, Flammer & Mendis (1986) using the KH theory. The 
KH instability theory has been successfully used to interpret the structure of 
the pc-scale jet in the radio source 3C273 (Lobanov & Zensus 2001). Birk et al. 
(2000) studied the role of KH instabilities in superwinds of primeval galaxies. 
The effect of radiative losses on the evolution of the KH instability in jets and 
outflows has been studied by many authors, cither in linear regime (e.g., Mas- 
saglia et al. 1992; Hardee & Stone 1997) or by direct numerical simulations 
(e.g., Rossi et al. 1997; Downes & Ray 1998). Recently, Michikoshi & Inutsuka 
(2006) studied the KH instability of the dust layer in protoplanetary discs to 
understand the effect of relative motion between gas and dust. 

However, in many astrophysical flows the single fluid approximation is in- 
appropriate. For example, the fact that the winds and the intergalactic medium 
may be only partially ionized should be considered in treatment of role of the 
KH instability in superwinds of primeval galaxies. We may have a similar situ- 
ation at the interface between a jet and its cocoon if the KH instability is to be 
viewed as a possible mechanism for entrainment. The problem becomes more 
complicated if we note that the development of the KH instability is strongly 
influenced by the magnetic field (e.g. Malagoli et al. 1996; Jones et al. 1997). 
Therefore it seems a two-fluid treatment of KH instability, taking account of 
the different motions of ions and neutrals and of the magnetic field, rather 
than a pure one fluid MHD description is a more appropriate approach. 

By neglecting the thermal pressure forces of the neutrals or perturbations 
in the neutral gas bulk velocities, the KH instability of a system consisting of 
ions and neutrals with applications to partially ionized cometary ionopauses, 
has been studied for both compressible (Ershkovich & Mendis 1986) and in- 
compressible cases (Ershkovich, Prialnik & Eviatar 1986; Ershkovich, Flammer 
& Mendis 1986). Chhajlani & Vyas (1991) incorporated the effects of finite 
resistivity and rotation on the KH instability of two superposed fluids slipping 
past each other. They showed that small rotation and the presence of neutral 
particles destabilize the system. Birk et al. (2000) studied the KH instability 
by taking into account the full incompressible dynamics of both the neutral 
and the ionized gas components with applications to multi-phase galactic out- 
flow winds. In another similar study, Watson et al. (2004) (hereafter WZHC) 
investigated the KH instability in the linear, partially ionized regime to de- 
termine its possible effect on entrainment in massive bipolar outflows. They 
showed that for much of the relevant parameter space, neutral and ions are 
sufficiently decoupled that the neutrals are unstable while ions are held in 
place by the magnetic field. Birk & Wicchcn (2002) focused on unstable shear 
flows in partially ionized dense dusty plasma. They considered the dust and 
neutral gas components so that dust and neutral collisions is the dominant mo- 
mentum transfer mechanism and dust component can interact with magnetic 
field lines, although dust charge fluctuations are negligible. They showed long 
wavelength modes can be stabilized by dust and neutral gas collisional mo- 
mentum transfer. By doing multifluid numerical simulations, Wiechen (2006) 
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studied KH modes in partially ionized, dusty plasmas for different masses and 
charges of dust. 

Analytical studies of KH instability in weakly ionized medium have so 
far considered one interface between two flows. For example, WHZC studied 
incompressible KH instability in a two-fluid system (charged and neutral) in 
order to understand entrainment in outflows from massive stars. They did 
the calculations in planar geometry with an interface at z = 0. However, it 
is not unreasonable to suppose that we may have a flowing layer with finite 
thickness and two interfaces. In this paper, we generalize the WHZC analysis 
by considering a finite thickness for a weakly ionised flow. We will show that 
the growth time-scale of the instability increases due to the thickness of the 
layer in particular in the long wavelength limit. In the next section, we will 
present the basic equations and assumptions of the model. In the section 3, 
we will discuss growing modes and their dependence on the thickness of the 
layer. 



2 General Formulation 

We follow the analysis of Chandrasekhar (1961) but take account of the dif- 
ferent bulk velocities and densities of the neutral and ions, both inside and 
outside the layer. These two components of our model are coupled by colli- 
sional interactions. We assume incompressibility for the analytical calculations 
and also that the convective term in the induction equation dominates the re- 
sistive one. For doing linear analysis, the unperturbed properties of the system 
are important. We suppose that the streaming takes place in the x— direction 
with velocity U(z), 

Tit \ f + U for 1*1 < d 
U W = \-U for|z|>d 

(U is constant) and the neutral and ion components have the same velocity in 
the unperturbed state. Note that 2d is the thickness of the layer and d > 0. The 
magnetic field is assumed parallel to the interface; that is, B = B x e x . Finally, 
all unperturbed physical quantities are assumed constant in each medium. The 
system under consideration consists of charged and neutral fluid components 
coupled by collisional interactions. The basic equations are 

V.v„,i - 0, (1) 

Pn(~Qf + ( V «- V ) v «) = -Vp„ - ■JniPnPitVn - V») (2) 

d v 1 

Pi(-fif + ( V *- V ) v = ^Pi + X B ) X B ~ lniPiPn{Vi - V n ) (3) 

ST) 

— = Vx(v 4 xB). (4) 

Also note that V.B = 0. In the above equations, 7„, is the collision rate 
coefficients per unit mass so that v ni — "f n iPi is the neutral-ion frequency. The 
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collision frequency determine the scale of the coupling between the different 
components and the coupling between each component and the magnetic field. 

Now, we perturb the physical variables as x{ z i x > t) = x'( z ) exp[i(u>t+k x x)], 
and so 

- k x p n u' n + iPn-r^ = 0, (5) 

<j)p n u' n = -k x p' n + ijniPnPi(u' n ~ u[), (6) 
<f>PnW' n = i^j- + HniPnPi(w' n - (7) 

where <fi = ui + k x U, and u' n and w' n are the x and the z components of the 
perturbed velocity of the neutrals. Also, the equations of the ions are 

- kxPiu'i + iPi-j 1 = 0, ( 8 ) 
<t>P%K = -k x Pi + i'Y ni p n p i (u' i - u' n ), (9) 

*»< = l ^-^^ + ^ri + * WiN ~ <>• (l0) 

where and w\ are the x and the z components of the perturbed velocity 
of the ions. We can reduce the above differential equations to a set of two 
differential equations for w\ and w' n as follows (WZHC) 

B 2 k 2 

[pi<t> ~ ~^f^ Vw 'i = HniPnPi'D(w' i - w' n ), (ll) 

p n 4>Vw' n = i'YniPnPiD{w' n - wl), (12) 

where V = 4h — k 2 . The behavior of the flow both inside and outside the layer 
is determined by the general solutions of the linear differential equations (|12p 
and (fTTj) . One can easily show that the general solutions of these equations 
are linear combinations of exp(+k x z) and exp(— k x z). We impose the condition 
that the perturbed quantities do not diverge at z — +00 and z = —00. Thus, 
the general solutions are 

C ie + k * z + C 2 e- k * z for < z < d 
^^H^-M i0Td < z ( 13 ) 



iC/'-HC^ 2 for0<z<d n .. 



C§e~ kxZ for d < z 

where Ci, C2, C3, C4, C5 and are constants which are determined by the 
boundary conditions. We note that, because of symmetry of the unperturbed 
state, there are solutions where Wi(z) = —Wi{—z) and w n (z) — —w n (—z) 
(odd solutions) and solutions where Wi(z) = Wi{—z) and w n (z) — w n (z) (even 
solutions). Thus, we can write the boundary condition at z = as 



Ci + sc 2 = 0, 



(15) 
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d + SC 5 = 0, (16) 

where S — +1 and S — — 1 are corresponding to even and odd solutions, re- 
spectively. Thus, considering symmetry properties of the solutions (i.e. equa- 
tions (fT5|) and (fl6|) ) our solutions are also valid for z < 0, if we transform 
z — ► — z. So, the solutions correctly vanish, as z — > — oo. 

Another boundary condition is the continuity of the ion and neutral dis- 
placements across the interfaces. So, the boundary conditions at z = d are 

de k * d + C 2 e~ k * d _ C 3 e- k * d 
u + k x U u> — k x U ' 

C 4 e k * d + C 5 e~ k * d _ C 6 e- k * d 
uj + k x ll ~ lu - k x U' 
Having equations (|15p. (|16p. (|17p and (fT5|) . we can rewrite Wi and io„ as 

„ x ( C 2 (e- k * z - Se k **) for0<z<d 
^W" \C 2 ^^(1-Se 2fe ^) ford<z (19) 

, ( s J C 5 (e- fc »* - Se fc * z ) for < z < d 

w n W-{ Cs (1 - 5e 2 ^ z ) for d < z (20) 

Clearly, the differential equations (fl2|) and (jlip are not valid at the in- 
terfaces z — d (or z = — d), where velocities of the species are not analytic. 
However, we can integrate each equation over an infinitesimal region around 
the interface z — d and use standard Gaussian pillbox arguments. Then we 
have i 

A = iA{ lniPn p^(w' n - wfi}, (21) 

A{ Pi 0—±} = iA{ lniPnPi — (wl - <)} + -^A{-— ±}, (22) 
dz dz 47r <b dz 



By substituting equations (fl9|) and (|20|) into the above equations and after 
lengthy (but straightforward) mathematical manipulations, we obtain 

+ m<e - (1 - X 2 )h 2 + a 2 ]C 2 - m£,C 5 = 0, (23) 

- £C 2 + {Z 2 + C - (1 - \ 2 )h 2 ]C 5 = 0, (24) 

where 

Pn , Uk x kvAi w 
m = — , ft = , a = , An — — 1£, (25) 

Pi V V V 

and A = S~ Y exp(— 2kd) and the ion Alfven velocity vj^i is defined as B/^Anpi. 
By the condition that the algebraic equations (f23|) and (|24|) have a nontrivial 
solution, we obtain the dispersion relation 

e + (m + l)e-[2(l-X 2 )h 2 -a 2 ]e 

- [(to + 1)(1 - X 2 )h 2 - a 2 }£ + (1 - A 2 )[(l - \ 2 )h 2 - a 2 ]h 2 = 0. (26) 
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This algebraic equation (|26|) is the basis of our stability analysis. As we see, 
the effect of the thickness of the layer appears by term A 2 which depends on 
the thickness of the layer. For a system with one interface at z = 0, WZHC 
obtained this algebraic equation 

£ 4 + (to + 1)£ 3 - (2h 2 - a 2 )£ 2 - [(to + l)h 2 - a 2 ]i + {h 2 - a 2 )h 2 = 0. 

We note that our algebraic equation |26|) reduces to the above equation if we 
set A = 0. 

3 Analysis 

Roots of equation (|26|) with positive real parts correspond to growing unstable 
modes. This equation is applicable for various weakly ionized systems. How- 
ever, in this work the input parameters are chosen to be appropriate for molec- 
ular outflows. This makes for easier comparison with the analysis of WZHC. In 
our illustrative examples, the relative velocity 2U is 20 Km s _1 , and the ratio, 
to, of neutral to ionized density is assumed to be around 10 6 . The collision fre- 
quency is v = 1.5 x 1CP 13 s _1 (Draine et al. 1983). The magnetic field strength 
is taken to be either 1 fj,G (weak-field case) or 1 mG (strong-field case). With 
these parameters the ratio of the Alfven velocity to the relative velocity is 
4400 and 4.4 in the strong and weak magnetic field cases, respectively. Having 
these input parameters, we can solve the algebraic equation (|26[) numerically 
to find the dispersion relation and the fastest growing modes. 

Figure Q] shows plots of the dispersion relation for the above input param- 
eters. When A — > 0, the dispersion relation tends to WZHC analysis (dashed 
curves in Figures[T]), as expected. Each curve is also labeled by the thickness of 
the layer, i.e. 2d. It can be seen from Figure[T]that the effect of the thickness of 
the layer is not significant for perturbations with short wavelength. However, 
as the wavelength of the perturbations increases, the thickness of the layer has 
a stabilizing effect on the growing modes and the dispersion relation deviates 
from the analysis of WZHC both in the weakly and strongly magnetised cases. 
Considering the thickness of the layer implies longer growth time-scale for the 
perturbations, in particular in the long wavelength limit. Moreover, as one 
would expect, as the thickness of the layer decreases the deviation from the 
WZHC analysis occurs at shorter wavelength, and also the growth time-scale 
becomes longer. 

Interestingly, the roots of equation (|26p which correspond to growing modes 
can be described by approximate analytical solutions. There are three approx- 
imate positive roots for this equation depending on its coefficients: 

(i) (1 - A 2 ) 1 ^ when (1 - X 2 )^ 2 h > 1; 

(ii) (1 - A 2 )/i 2 , when (1 - X 2 )^ 2 h < 1 and (1 - A 2 ) 1 / 2 < {a/h)m,- 1 / 2 - 

(iii) (1 - X 2 f/ 2 h(l - {a/h) 2 (l - A 2 )- 1 /™" 1 ) 172 and (1 - A 2 ) 1 / 2 > {a/h^ 1 ' 2 . 
The agreement between these approximate roots and numerical roots of the 

equation (|26|) is excellent. Thus, as long as (1 — A 2 ) 1 / 2 /! > 1 the only growing 
mode is the first growing mode, irrespective of the strength of the magnetic 
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Fig. 1 Growth rate of the instability vs. wavelength of the perturbations for two cases 
B = 1 fj,G (top) and B = 1 mG (bottom). In both cases, it is assumed U = 10 Km s -1 . 
The dashed curves correspond to the WZHC analysis and the other curves are labeled by 
the thickness of layer. Black circles in the bottom plot denote the wavelength at which the 
growing mode changes from (1 — X 2 ) 1 ^ 2 h to (1 — X 2 )h 2 . 
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field (see top plot of Figure [I]). Once this inequality is violated, we may have 
the second growing mode, depending on whether (1 — A 2 ) 1 / 2 < {a/lijmr 1 / 2 . 
For example, for B = 1 mG we have (a//i)m -1 / 2 = 4.5 > 1 > (1 — A 2 ) 1 / 2 
and so the transition from the first to the second growing modes occurs when 
(1 — \ 2 ) l / 2 h = 1. We can denote this transition wavelength by At and is shown 
by black circles in Figure [T] (bottom plot). But in the weakly magnetic case 
where B = 1 [iG, we have (a//i)m- 1 / 2 = 4.5 x 1CT 3 < (1 - A 2 ) 1 / 2 . So, for 
wavelengths greater than A*, the third growing mode appears, although the 
second term inside the parenthesis of this root is much smaller than unity and 
the third mode actually tends toward first mode (top plot). Considering our 
ranges of the input parameters, the third mode is very close to the first mode. 

The first mode corresponds to the usual hydrodynamic instability, but the 
second mode appears in our weakly ionized case as a new mode (WZHC). The 
transition between these two modes depends on the thickness of the layer. 
Figure [2] shows the transition wavelength At versus the the half thickness of 
the layer. For thin layers, the transition wavelength is smaller than the thick 
layers. In fact, as the thickness of the layer increases, this transition wavelength 
increase and tends to an asymptotic value for the very thick layers. WZHC 
analysis corresponds to this asymptotic regime. 



4 Conclusions 

In this study we have extended the analysis of WZHC by considering a layer of 
ions and neutrals with finite thickness interacting with an ambient medium in 
the incompressible limit. However, this model can not be applied directly for 
systems, such as outflows or jets because of our simplifying assumptions. We 
can, however, gain insight into the possible effects of shear between a moving 
layer of ions and neutrals and the surrounding medium. 

To our knowledge, there are a few numerical simulations of KH instability 
of a multifluid system (Birk et al. 2000; Birk & Wiechen 2002; Wiechen 2006). 
While Birk et al. (2000) showed that KH modes can operate fast enough to 
amplify the magnetic field strengths in superwinds of primeval galaxies within 
the timescale of the outflow, Birk & Wiechen (2002) and Wiechen (2006) 
studied KH instability in a dusty magnetized medium. It is difficult to make a 
direct comparison between our results and these numerical simulations because 
of different basic assumptions. We considered an incompressible two-fluid layer, 
but all these numerical simulations are for compressible systems. Moreover, the 
effect of charged dust particles are considered in Birk & Wiechen (2002) and 
Wiechen (2006), though we have neglected charged grains in our analysis. 

We have shown that the growth time-scale of KH instability in a layer of 
ions and neutrals increases when compared to a system with one interface, 
this effect being more evident in the long wavelength limit. There are two 
dominant growing modes of KH instability in such a system, depending on 
the magnetic strength and the thickness of the layer. As the thickness of the 
layer decreases, the transition between these two unstable modes occurs at 
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Fig. 2 The transition wavelength At vs. the half thickness d of the layer. 



shorter wavelengths. Given these results, it remains to future work to extend 
this analysis to compressible case and make direct comparisons with jets or 
outflows with finite thickness. 
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